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ABSTRACT 

Magnetohydrodynamic (MHD) turbulence is encountered in a wide variety of astrophysical plasmas, including 
accretion disks, the solar wind, and the interstellar and intracluster medium. On small scales, this turbulence is 
often expected to consist of highly anisotropic fluctuations with frequencies small compared to the ion cyclotron 
frequency. For a number of applications, the small scales are also collisionless, so a kinetic treatment of the 
turbulence is necessary. We show that this anisotropic turbulence is well described by a low frequency expansion 
of the kinetic theory called gyrokinetics. This paper is the first in a series to examine turbulent astrophysical 
plasmas in the gyrokinetic limit. We derive and explain the nonlinear gyrokinetic equations and explore the linear 
properties of gyrokinetics as a prelude to nonlinear simulations. The linear dispersion relation for gyrokinetics is 
obtained and its solutions are compared to those of hot-plasma kinetic theory. These results are used to validate the 
performance of the gyrokinetic simulation code GS2 in the parameter regimes relevant for astrophysical plasmas. 
New results on global energy conservation in gyrokinetics are also derived. We briefly outline several of the 
problems to be addressed by future nonlinear simulations, including particle heating by turbulence in hot accretion 
flows and in the solar wind, the magnetic and electric field power spectra in the solar wind, and the origin of small- 
scale density fluctuations in the interstellar medium. 

Subject headings: Magnetic fields — methods: analytical — methods: numerical — plasmas — turbulence 



1. INTRODUCTION 

The Goldreich and Sridhai- (1995) (hereafter GS) theory of 
MHD turbulence (see also Sridhai" and Goldreich 1994; Gol- 
dreich and Sridhar 1997; Li th wick and Goldreich 2001) in a 
mean magnetic field predicts that the energy cascades primar- 
ily by developing small scales perpendicular to the local field, 
with k± ^ as schematically shown in Figure 1 (cf. earlier 
work by Montgomery and Turner 1981; Shebalin et al. 1983). 
Numerical simulations of magnetized turbulence with a dynam- 
ically strong mean field support the idea that such a turbulence 
is strongly anisotropic (Maron and Goldreich 2001; Cho et al. 

2002) . In situ measurements of turbulence in the solar wind 
(Belcher and Davis 1971; Matthaeus et al. 1990) and observa- 
tions of interstellar scintillation (Wilkinson et al. 1994; Trotter 
et al. 1998; Rickett et al. 2002; Dennett-Thorpe and de Bruyn 

2003) also provide evidence for significant anisotropy. 

In many astrophysical environments, small-scale perturba- 
tions in the MHD cascade have (parallel) wavelengths much 
smaller than, or at least comparable to, the ion mean free path, 
implying that the turbulent dynamics should be calculated using 
kinetic theory. As a result of the intrinsic anisotropy of MHD 
turbulence, the small-scale perturbations also have frequencies 
well below the ion cyclotron frequency, oj ^ fi,, even when the 
perpendicular wavelengths are comparable to the ion gyrora- 
dius (see Figure 1). Anisotropic MHD turbulence in plasmas 
with weak collisionality can be described by a system of equa- 
tions called gyrokinetics. 

Particle motion in the small-scale turbulence is dominantly 



the cyclotron motion about the unperturbed field lines. Gyroki- 
netics exploits the time-scale separation (w ^ il,) for the elec- 
tromagnetic fluctuations to eliminate one degree of freedom in 
the kinetic description, thereby reducing the problem from 6 to 
5 dimensions (three spatial plus two in velocity space). It does 
so by averaging over the fast cyclotron motion of charged par- 
ticles in the mean magnetic field. The resulting "gyrokinetic" 
equations describe charged "rings" moving in the ring-averaged 
electromagnetic fields. The removal of one dimension of phase 
space and the fast cyclotron time-scale achieves a more com- 
putationally efficient treatment of low frequency turbulent dy- 
namics. The gyrokinetic system is completed by electromag- 
netic field equations that are obtained by applying the gyroki- 
netic approximation to Maxwell's equations. The gyrokinetic 
approximation orders out the fast MHD wave and the cyclotron 
resonance, but retains finite-Larmor-radius effects, collisionless 
dissipation via the parallel Landau resonance, and collisions. 
Both the slow MHD wave and the Alfven wave are captured by 
the gyrokinetics, though the former is damped when its wave- 
length along the magnetic field is smaller than the ion mean free 
path (Barnes 1966; Foote and Kulsrud 1979). 

Linear (Rutherford and Frieman 1968; Taylor and Hastie 
1968; Catto 1978; Antonsen and Lane 1980; Catto et al. 1981) 
and nonlinear gyrokinetic theory (Frieman and Chen 1982; Du- 
bin et al. 1983; Lee 1983, 1987; Hahm et al. 1988; Brizard 
1992) has proven to be a valuable tool in the study of labo- 
ratory plasmas. It has been extensively employed to study the 
development of turbulence driven by microinstabilities, e.g., the 
ion- and electron-temperature-gradient instabilities (e.g., Dim- 
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its et al. 1996; Dorland et al. 2000; Jenko et al. 2000; Rogers 
et al. 2000; Jenko et al. 2001; Jenko and Dorland 2001, 2002; 
Candy et al. 2004; Parker et al. 2004). For these applications, 
the structure of the mean equilibrium field and the gradients in 
mean temperature and density are critical. The full gyrokinetic 
equations allow for a spatially varying mean magnetic field, 
temperature, and density. In an astrophysical plasma, the mi- 
croinstabilities associated with the mean spatial gradients are 
unlikely to be as important as the MHD turbulence produced 
by violent large scale events or instabilities. Our goal is to use 
gyrokinetics to describe this MHD turbulence on small scales 
(see Fig. 1). For this purpose, the variation of the large-scale 
field, temperature, and density is unimportant. We, therefore, 
consider gyrokinetics in a uniform equilibrium field with no 
mean temperature or density gradients. 
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Fig. 1. — Schematic diagram of the low-frequency, anisotropic Alfv6n-wave 
cascade in wavenumber space: the horizontal axis is perpendicular wavenum- 
ber; the vertical axis is the parallel wavenumber, proportional to the frequency. 
MHD (or, rather, strictly speaking, its Alfv6n-wave part; see Schekochihin et 
al. 2006) is valid only in the limit w <S Q, and k±pi <C 1; gyrokinetic theory 
remains vaUd when the perpendicular wavenumber is of the order of the ion 
Larmor radius, k±pi ~ 1. Note that u) — > H, only when fcyp,- — » 1, so gyroki- 
netics is appUcable for fcy <g k±. 

This is the first in a series of papers to apply gyrokinetic the- 
ory to the study of turbulent astrophysical plasmas. In this pa- 
per, we derive the equations of gyrokinetics in a uniform equi- 
librium field and explain their physical meaning. We also derive 
and analyze the linear gyrokinetic dispersion relation in the col- 
lisionless regime, including the high-beta regime, which is of 
particular interest in astrophysics. Future papers will present 
analytical reductions of the gyrokinetic equations in various 
asymptotic limits (Schekochihin et al. 2006) and nonlinear sim- 
ulations applied to specific astrophysical problems These sim- 
ulations use the gyrokinetic simulation code GS2. As part of 
our continuing tests of GS2, we demonstrate here that it ac- 
curately reproduces the frequencies and damping rates of the 
linear modes. 

The remainder of this paper is organized as follows. In §2, 



we present the gyrokinetic system of equations, giving a phys- 
ical interpretation of the various terms and a detailed discus- 
sion of the gyrokinetic approximation. A detailed derivation 
of the gyrokinetic equations in the limit of a uniform, mean 
magnetic field with no mean temperature or density gradients 
is presented in Appendix A. Appendix B contains the first pub- 
lished derivation of the long-term transport and heating equa- 
tions that describe the evolution of the equiUbrium plasma — 
these are summarized in §2.5. In §2.6, we introduce the linear 
collisionless dispersion relation of the gyrokinetic theory (de- 
tailed derivation given in Appendix C, various analytical hmits 
worked out in Appendix D). The numerics are presented in 
§3, where we compare the gyrokinetic dispersion relation and 
its analytically tractable limits with the full dispersion relation 
of the hot-plasma theory and with numerical results from the 
gyrokinetic simulation code GS2. We also discuss the effect 
of colUsions on collisionless damping rates (§3.5) and illus- 
trate the limits of applicability of the gyrokinetic approximation 
(§3.6). Finally, §4 summarizes our results and discusses several 
potential astrophysical applications of gyrokinetics. Definitions 
adopted for plasma parameters in this paper are given in Ta- 
ble 1. 

2. GYROKINETICS 

In this section we describe the gyrokinetic approximation and 
present the gyrokinetic equations themselves in a simple and 
physically motivated manner (the details of the derivations are 
given in the Appendices). To avoid obscuring the physics with 
the complexity of gyrokinetics in full generahty, we treat the 
simplified case of a plasma in a straight, uniform mean mag- 
netic field. Bo = BqZ and with a spatially uniform equilibrium 
distribution function, VFq = (the slab limit). This case is also 
of most direct astrophysical relevance because the mean gradi- 
ents in turbulent astrophysical plasmas are generally dynami- 
cally unimportant on length scales comparable to the ion gyro- 
radius. 

2.1. The Gyrokinetic Ordering 

The most basic assumptions that must be satisfied for the gy- 
rokinetic equations to be apphcable are weak couphng, strong 
magnetization, low frequencies and small fluctuations. The 
weak coupling is the standard assumption of plasma physics: 
noe^oe ^ 1' where riQe is the mean electron number density and 
Xoe is the electron Debye length. This approximation allows 
the use of the Fokker-Planck equation to describe the kinetic 
evolution of all plasma species. 

The conditions of strong magnetization and low frequencies 
in gyrokinetics mean that the ion Larmor radius must be 
much smaller than the macroscopic length scale L of the equi- 
librium plasma and that the frequency of fluctuations ui must be 
small compared to the ion cyclotron frequency 12,, 

Pi=^^L, (1) 

The latter assumption allows one to average all quantities over 
the Larmor orbits of particles, one of the key simplifications al- 
lowed by the gyrokinetic theory. Note that the assumption of 
strong magnetization does not require the plasma beta (the ra- 
tio of the thermal to the magnetic pressure, f3 = Enp/B^) to be 
small. A high-beta plasma can satisfy this constraint as long as 
the ion Larmor radius is small compared to the gradients of the 
equilibrium system. In most astrophysical contexts, even a very 
weak magnetic field meets this requirement. 
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Table 1 

Definitions of plasma parameters. 







m 


nnrtiplp ma*;*; 


Hs 


nartiplp fhartrp iw/f tak'p n- — — n — /^^ 


"'VS 


number density 




temperature (in units of energy) 


= 8™q.sTo, / Bl 


plasma beta 


Vrt, =(2ros/m,)'/2 


thermal speed 


c, = (ro,/m,)'/2 


sound speed 


VA = Bo/(47rm,no;)'/^ 


Alfven speed 


c 


speed of light 


= q,Bo/(m,c) 


cyclotron frequency (carries sign of q,) 




Lamior radius 



To derive the gyrokinetic equations, we order the time and 
length scales in the problem to separate fluctuating and equi- 
librium quantities. The remainder of this section defines this 
formal ordering and describes some simple consequences that 
follow from it. 

Two length scales are characteristic of gyrokinetics: the 
small length scale, the ion Larmor radius p,, and the larger 
length scale Iq, which is here introduced formally and will be 
argued below to be the typical parallel wavelength of the fluctu- 
ations. Their ratio defines the fundamental expansion parameter 
e used in the formal ordering: 

Pi 



«: 1. 



(2) 



There are three relevant time scales, or frequencies, of in- 
terest. The fast time scale is given by the ion cyclotron fre- 
quency il, . The distribution function and the electric and mag- 
netic fields are assumed to be stationary on this time scale. The 
intermediate time scale corresponds to the frequency of the tur- 
bulent fluctuations. 



^ ~ 0(en,)- 
to 



(3) 



The slow time scale is connected to the rate of heating in the 
system, ordered as foUows 

_L^e2^^~0(e3f2,). (4) 

theat 'o 

The distribution function / of each species s (= e,;, the 
species index is omitted unless necessary) and magnetic and 
electric fields B and E are split into equilibrium parts (denoted 
with a subscript 0) that vary at the slow heating rate and fluctu- 
ating parts (denoted with 5 and a subscript indicating the order 
in e) that vary at the intermediate frequency lo: 

f(r,\j) = Fn(\, t) + (5/i(r , V, f ) + (5/2(r, v, f) + • • • , (5) 
B(r,0 = Bo + 5B(r,f) = Boz + Vx A, (6) 
1 (9A 

E(r,O = <5E(r,f) = -V0-- — . (7) 
c at 

Let us now list the gyrokinetic ordering assumptions. 

• Small fluctuations about the equilibrium. Fluctuating 
quantities are formaUy of order e in the gyrokinetic ex- 
pansion, 

5fi (5B (5E 



0{e). 



(8) 



Bq (vthJc)BQ 
Note that although fluctuations are small, the theory is 
fully nonlinear (interactions are strong). 



• Slow-time-scale variation of the equilibrium. The equi- 
librium varies on the heating time-scale. 



]_dFo 

Fq dt 



O 



1 

^hea 



■0{e 



(9) 



Derivations for laboratory plasmas (Frieman and Chen 
1982) have included a large-scale [~ 0(1 /Iq)] spatial 
variation of the equilibrium (Fq and Bo) — this we omit. 
The slow-time-scale evolution of the equilibrium, how- 
ever, is treated for the first time here. 



Intermediate-time-scale variation of the fluctuating 
quantities. The fluctuating quantities vary on the inter- 
mediate time scale 



1 dSf 1 ajB 1 ajE fvth, 

'^^jf~dr^W\~dr^\5K\~w^ \k 



(10) 



Intermediate-time-scale collisions. The coUision rate in 
gyrokinetics is ordered to be the same as the intermedi- 
ate time-scale 



zy~o|^y^ ) ~0(a;). 



(11) 



Colhsionless dynamics with a; > z/ are treated correctly 
as long as 1/ > ew. 

• Small-scale spatial variation of fluctuations across the 
mean fleld. Across the mean magnetic field, the fluctua- 
tions occur on the small length scale 



z X V(5/ z X V(5B z X V5E ^ / 1 



Sf 



\6B\ 



Pi 



(12) 



• Large-scale spatial variation of fluctuations along the 
mean fleld. Along the mean magnetic field the fluctua- 
tions occur on the larger length scale 



z-VSf z-VJB Z-V6E 



Sf 



\5B\ 



0{ 



(13) 
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Fig. 2. — Diagram illustrating the ring average in gyrokinetics. The ion 
position is given by the open circle, the guiding center position by the filled 
circle; the ring average, centered on the guiding center position, is denoted by 
the thick-lined circle passing through the particle's position. The characteristic 
perpendicular and parallel length scales in gyrokinetics are marked and /|| ; 
here the perpendicular scale is exaggerated for clarity. The unperturbed mag- 
netic field Bq is given by the long-dashed line and the perturbed magnetic field 
B by the solid line. The particle drifts off of the field line at ux, roughly the 
E X B velocity. 



With a small length scale across the field and a large length 
scale along the field, the typical gyrokinetic fluctuation is 
highly anisotropic: 

h 



^ ~ ^ ~ 0(e). 



(14) 



Figure 2 presents a schematic diagram that depicts the length- 
scales associated with the gyrokinetic ordering. The typical per- 
pendicular flow velocity, roughly the E x B velocity, is 

c^ExB , 
u_L ~ — Oievth)- (15) 

^0 

The typical perpendicular fluid displacement is ~ m_l/w ~ 
1/^_L ~ 0{pi), as is the field line displacement. Since displace- 
ments are of order the perpendicular wavelength or eddy size, 
the fluctuations are fully nonlinear. 

2.2. The Gyrokinetic Ordering andMHD Turbulence 

The GS theory of incompressible MHD turbulence conjec- 
tures that, on sufficiently small scales, fluctuations at all spa- 
tial scales always arrange themselves in such a way that the 
Alfven time scale and the nonlinear decorrelation time scale 
are similar, uj k\\VA ^ k±_u_\_. This is known as the critical 
balance. A modification of Kolmogorov (1941) dimensional 
theory based on this additional assumption then leads to the 
scahng m_l ~ U{k±_L)~^l^, where U and L are the velocity and 
the scale at which the turbulence is driven. For a detailed dis- 
cussion of these results, we refer the reader to Goldreich and 
Sridhar's original papers or to a review by Schekochihin and 
Cowley (2006a). 

Here, we show that the gyrokinetic ordering is manifestly 
consistent with, and indeed can be constructed on the basis of. 



the GS critical balance conjecture. Using the critical balance 
and the GS scahng of m^, we find that the ratio of the turbulent 
frequency w ^k^u^ to the ion cyclotron frequency is 

^-^(W/^(f)^^ (16) 
The ratio of parallel and perpendicular wave numbers is 

&~^»->-"(£)"- 

Both of these ratios have to be order e in the gyrokinetic expan- 
sion. Therefore, for the GS model of magnetized turbulence, 
we define flie expansion parameter 

-(!)"'• 

Comparing fliis with equation (2), we may formally define the 
length scale Zq used in the gyrokinetic ordering as Iq = p'J^V-l^ . 
Physically, this definition means that /q is the characteristic par- 
allel length scale of the turbulent fluctuations in the context of 
GS turbulence. Note that our assumption of no spatial variation 
of the equilibrium is, therefore, equivalent to assuming that the 
variation scale of Fq and Bq is Iq — ^this is satisfied, e.g., for 
the injection scale L. 

One might worry that the power of 1/3 in equation (18) 
means that the expansion is valid only in extreme circum- 
stances. For astrophysical plasmas, however, is so small 
that this is probably not a significant restriction. To take the 
interstellar medium as a concrete example, p, ~ 10^ cm and 

1 /% 

L ~ 100 pc ~ 3 X 10^° cm (the supernova scale), so (pi/L) 

1 /3 

10"^. For galaxy clusters, (pi/L) ~ 10"^ (Schekochihin and 
Cowley 2006b), for hot, radiatively inefficient accretion flows 

around black holes, ^ 10~^ (Quataert 1998), and for 

1 /3 

the solar wind, (pi/L) ~ 10"^. 

In the gyrokinetic ordering, all cyclotron frequency effects 
(such as the cyclotron resonances) and the fast MHD wave are 
ordered out (for a more general approach using a kinetic de- 
scription of plasmas in the gyrocenter coordinates that retains 
the high-frequency physics, see the gyrocenter-gauge theory of 
Qin et al. 2000). The slow and Alfven waves are retained, 
however, and collisionless dissipation of fluctuations occurs via 
the Landau resonance, through Landau damping and transit- 
time, or Barnes (1966), damping. The slow and Alfven waves 
are accurately described for arbitrary k±pi, as long as they are 
anisotropic ik\\ ~ ek±). Subsidiary ordering of collisions (as 
long as it does not interfere with the primary gyrokinetic or- 
dering) allows for a treatment of collisionless and/or coUisional 
dynamics. Similarly, subsidiary ordering of the plasma beta al- 
lows for both low- and high-beta plasmas. 

The validity of GS turbulence theory for compressible astro- 
physical plasmas is an important question. Direct numerical 
simulations of compressible MHD turbulence (Cho & Lazar- 
ian 2003) demonstrate that spectrum and anisotropy of slow 
and Alfven waves are consistent with the GS predictions. A 
recent work exploring weak compressible MHD turbulence in 
low-beta plasmas (Chandran 2005) shows that interactions of 
Alfven waves with fast waves produce only a small amount of 
energy at high k\\ (weak turbulence theory for incompressible 
MHD predicts no energy at high fcy), but this is unlikely to alter 
the prediction of anisotropy in strong turbulence. Both of these 
works demonstrate that a small amount of fast- wave energy cas- 
cades to high frequencies, but the dynamics in this regime are 
energetically dominated by the low-frequency Alfven waves. 



ASTROPHYSICAL GYROKINETICS 



5 



2.3 . Coordinates and Ring Averages 

Gyrokinetics is most naturally described in guiding center 
coordinates, where the position of a particle r and the position 
of its guiding center are related by 

V X z 



r = R. 



(19) 



The particle velocity can be decomposed in terms of the par- 
allel velocity vy, perpendicular velocity v^, and the gyrophase 
angle 6: 

V = V|[Z + vj^(cos6'x + sin6'y). (20) 

Gyrokinetics averages over the Larmor motion of particles 
and describes the evolution of a distribution of rings rather than 
individual particles. The formalism requires defining two types 
of ring averages: the ring average at fixed guiding center R,, 



(fl(r,v,f))R, = 



dOaiR, 



V X z 



-,v,f), 



(21) 



where the 9 integration is done keeping R^, v± and V|| constant, 
and the ring average at fixed position r, 

1 

2^ 



(fl(R.,V,f))r 



V X z 

d9a(r+— — ,y,t), 



(22) 



where the integration is at constant r, v_l, vy . 

2.4. The Gyrokinetic Equations 

The detailed derivation of the gyrokinetic equations is given 
in Appendix A. Here we summarize the results and their phys- 
ical interpretation. The full plasma distribution function is ex- 
panded as follows 

qs(l>{r,ty 



fs = Fosiv,t)exp 



Tos 



+hsiRs,v,v±,t) + Sf2s + --- , (23) 



where v = (v^ + v|)^/^ and the equilibrium distribution function 
is Maxwellian: 

,.2 \ 



exp I -— 



v; 



(24) 



The first-order part of the distribution function is com- 
posed of a term that comes from the Boltzmann factor, 
exp[-^.5(^(r,f)/roj| ^ \-qs(f)(r,t)/T{)s, and the ring distribution 
hs- The ring distribution /z, is a function of the guiding center 
position Rj (not the particle position r) and two velocity coor- 
dinates, V and vx-^ It satisfies the gyrokinetic equation: 



dh, , dh, c r 



dlh 
'dt 



coU 



<9(x)r. Fa, 



dt 



Tqs 

(25) 

where the electromagnetic field enters via the ring average of 
the gyrokenetic potential x = (/) - v • A/c. The Poisson bracket 
is defined by [U,V] =z-[{dU/dRs) x (dV/dRs)]. The scalar 
potential ^ and the vector potential A are expressed in terms 
of hs via Maxwell's equations: the Poisson's equation, which 
takes the form of the quasineutrality condition 



E 



(t> + q, / d\{h,)r]=0. 



(26) 



the parallel component of Ampere's law, 

s ^ 

^ Note that in the inhomogeneous case, it is more convenient to use the energy m, 
instead of v and vx- 



and the perpendicular component of Ampere's law, 

Vx(5fi|| = j dMi^ X y±)hs)r- (28) 

The gyrokinetic equation (25) can be written in the follow- 
ing, perhaps more physically illuminating, form. 





dRs\ 


dh, 


(dh, 


dt \ 


dt 1 


K,' dRs 


V dt 


where 








/dR, 




c 




\ dt 


)r. 




dRs 



= VllZ- 



c d{(l))R 
Bo dR, 



coU 



X Z 



X z 



d£s 



Fos 



dt I Tqs 



(29) 



is the ring velocity. Eg = /'2)msV^ + qs<p is the total energy of 
the particle and 

'd£A _ 5(x)r, 
dt ' "^^ 



dt 



(31) 



(32) 



Note that the right-hand side of equation (29) is 
Fo^^_/d£^dfA 
dt I R, \ dt d£s I ^ 

written to lowest order in e. Using this and the conserva- 
tion of the first adiabatic invariant, {dfis/dt)^^ = 0, where fig = 
msv\/2BQ, it becomes clear that equation (29) is simply the 
gyroaveraged Fokker-Planck equation 

'f-(t)J.- - 

where only the lowest order in e has been retained. 

A simple physical interpretation can now be given for each 
term in equation (29). It describes the evolution of a distribution 
of rings hs that is subject to a number of physical influences: 

• motion of the ring along the ring-averaged total (per- 
turbed) magnetic field: since VA|| x z = (5Bx, 



dRs Bo J dRs 



B 



dhs 
dRs 



the ring averaged E x B drift: 

C 5(<^)r, 



Bo dRs 

the VB drift: 
'(j (vx ■ Ax)r, 
dRs 



dhs _/ E X Bo 

dR:~v~Br 



dhs 



(34) 



; (35) 



fin 



dhs 
dR. 



dh. 



(5B|| 

Bo / T, dRs ' 



(36) 

where, if we expand the ring average [equation (21)] in 
small v X z/fls, we get, to lowest order, the familiar drift 
velocity: 

'SBi 



I' SB X Bp 
c — , (37) 



Bo '^/^^ ' qsBl 

where /Zj = msv\/2Bo is the first adiabtic invariant 
(magnetic moment of the ring) and V5 = V^fiy is taken 
at the center of the ring: r = R,; 

and the first adiabatic invariant (the magnetic moment) iXs = (l/2)msVj_/Bo 
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• the (linearized) effect of collisions on the perturbed ring 
distribution function: —{dhs/ 9f)coii (the gyrokinetic col- 
lision operator is discussed in detail in Schekochihin et 
al. 2006); 

• the effect of coUisionless work done on the rings by the 
fields (the wave-ring interaction): the right-hand side of 
equation (29). 

We have referred to the ring averaged versions of the more 
familiar guiding center drifts. Figure 2 shows the drift of the 
ring along and across the magnetic field. 

2.5. Heating in Gyrokinetics 

The set of equations given in the previous section determines 
the evolution of the perturbed ring distribution and the field 
fluctuations on the intermediate time-scale characteristic of the 
turbulent fluctuations. To obtain the evolution of the distribu- 
tion function Fq on the slow (heating) time-scale, we must con- 
tinue the expansion to order e^. This is done in Appendix B, 
where the derivations of the particle transport and heating equa- 
tions for our homogeneous equilibrium, including the equation 
defining the conservation of energy in externally driven systems 
(e.g., "forced" turbulence), are given for the first time. In an 
inhomogeneous plasma, turbulent diffusion, or transport, also 
enters at this order and proceeds on the slow time scale. Let us 
summarize the main results on heating. 

In the homogeneous case, there is no particle transport on the 
slow time-scale, 

dnps 
dt 

The evolution of the temperature Tqj of species s on this time 
scale is given by the heating equation 



resolve the entropy production. The numerical demonstration 
of the collisional heating and its independence of the collision 
rate is given in §3.5. 

In the homogeneous case, turbulence will damp away un- 
less driven. In our simulations, we study the steady-state ho- 
mogeneous turbulence driven via an external antenna current 
ja introduced into Ampere's law — i.e., the parallel and perpen- 
dicular components of ja are added to the right-hand sides of 
equations (27)-(28). The work done by the antenna satisfies 
the power-balance equation 



d \ — 



Fos 



\ dt 



hs{^] ) (40) 



coll 



= 0. 



(38) 



(see Appendix B.3). Thus, the energy input from the driving 
antenna is dissipated by heating the plasma species. The lesson 
of equation (39) is that this heat is always produced by entropy- 
increasing colhsions. 

2.6. Linear CoUisionless Dispersion Relation 

The derivation of the linear dispersion relation from the gy- 
rokinetic equations (25)-(28) is a straightfoward hnearization 
procedure. In Appendix C, it is carried out step by step. A key 
technical fact in this derivation is that once the electromagnetic 
fields and the gyrokinetic distribution function are expanded in 
plane waves, the ring averages appearing in the equations can 
be written as multiplications by Bessel functions. The resulting 
dispersion relation for Unear, coUisionless gyrokinetics can be 
written in the foUowing form 



aj2 



-AB+B' 



= iAE+BQ'- 



2%. 



dip, 
dt 



dS 



V 



-qs 



dt 



hs + nosi^'EiTor-Tos) 



dh 



,3 Tqj 



h. 



dfh 
'dt 



coll 



where w = io/\k\\ \va and, taking qi = -qe = e, not = noe, 

he 



B=\ 



(39) 



-ro(a,)+^[i- 



ro(a<,)], 



(41) 



(42) 



(43) 



D = 2ri(a,)eZ(e) + 2^ri(a,)6Z(6), 
£ = ri(a,)-ri(a,), 



(44) 



(45) 



(46) 



The overbar denotes the medium-time average over time Af 
such that l/w <C <C 1/ (e^w) [see equation (B 1)]. The second 
term on the right-hand side (proportional to i^J) corresponds to 
the collisional energy exchange (see, e.g., Helander & Sigmar 
2002) between species r and s^. It is clear from the second 
equality in equation (39) that the heating is ultimately always 
colhsional, as it must be, because entropy can only increase 
due to collisions. When the coUisionality is small, u <^ui,\he 
heating is due to the coUisionless Landau damping in the sense 
that the distribution function develops small-scale structure 
in velocity space, with velocity scales Av ^ (9(i^'''^). Colh- 
sions smooth these small scales at the rate vvji^. / (Av)^ ~ w, so 
that the heating rate [given by the second expression in equa- 
tion (39)] becomes asymptotically independent of jy in the coUi- 
sionless hmit (see related discussion of colhsionless dissipation 
by Krommes and Hu 1994; Krommes 1999). We stress that it is 
essential for any kinetic code, such as GS2, to have some col- 
lisions to smooth the velocity distributions at small scales and 

* We have been cavalier about treating the collision operator up to this point. The characteristic time scale of the interspecies collisional heat exchange is 
~ i''\me/mi)^l^(Jt)i-TQe)lToe- For the two terms on the right-hand side of equation (39) to be formally of the same order, we must stipulate i/"(mc/mj)'^^(7bi- 
1'oe)/Toe ~ 0(e^u). This ordering not only ensures that the zeroth-order distribution function is a Maxwellian but also provides greater flexibility in ordering the 
coUisionality relative to the intermediate time-scale of the fluctuations. We ignore this technical detail here — in most cases, the second term on the right-hand side of 
equation (39) is small compared to the first term, allowing a relatively large temperature difference between species to be maintained. 



where = co/\k\\\vth,, Z{^s) is flie plasma dispersion func- 
tion, a, = k\p^/2, and ro(a,v) = I()(as)e~"% ri(Q;5) = [loia,)- 
/i(Q;j)]e~"' (/q and /i are modified Bessel functions). These 
functions arise from velocity-space integrations and ring aver- 
ages: see Appendix C for details. 

The complex eigenvalue solution U to the equation (41) de- 
pends on three dimensionless parameters: the ratio of the ion 
Larmor radius to the perpendicular wavelength, k±pi; the ion 
plasma beta, or the ratio of ion thermal pressure to magnetic 
pressure, /?,; and the ion to electron temperature ratio, Tbi/^e- 
Thus, uJ = ojGK(k±Pi,/3i, Toi/Toe). 
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2.6. 1 . Long-Wavelength Limit 

Let us first consider the linear physics at scales large com- 
pared to the ion Larmor radius, where the comparison to MHD 
is more straightforward. These are not new results, but they 
are an important starting point for the more general results to 
follow. First, recall the MHD waves in the anisotropic Umit 



to : 



± 



Alfven waves, 
slow waves, 



(47) 
(48) 







fast magnetosonic waves, (49) 

entropy mode, (50) 

where c, is the sound speed. The fast magnetosonic waves have 
been ordered out of gyrokinetics because, when k_\_pi ~ 1, their 
frequency is of order the cyclotron frequency f2,. The removal 
of the fast waves is achieved by balancing the perpendicular 
plasma pressure with the magnetic field pressure [see equa- 
tion (A30)]. Here we are concerned with the Alfven and slow 
waves in the colhsionless hmit. Note that the entropy mode 
is mixed with the slow-wave mode when the parallel wave- 
length is below the ion mean free path. In this paper, when- 
ever we refer to the "slow-wave" part of the dispersion relation, 
we sacrifice terminological precision to brevity. Strictly speak- 
ing, the "slow waves," as understood below, are everything that 
is not Alfven waves — namely, modes involving fluctuations of 
the magnetic-field strength, which can also be aperiodic (have 
zero real frequency) (for further discussion of this component 
of gyrokinetic turbulence, see Schekochihin et al. 2006). 

The left-hand side of equation (41) contains two factors. We 
will see that the first factor corresponds to the Alfven-wave 
solution, the second to the slow-wave solution. The right- 
hand side of equation (41) represents the coupling between the 
Alfven and slow waves that is only important at finite ion Lar- 
mor radius. 

In the long- wavelength limit k±pi ^ 1, or a,- ^ 1, we can 
expand ro(as) — 1 -a,s and ri(Qj) ~ 1 -3q:.5/2. We can also ne- 
glect terms that multiply powers of the electron-ion mass ratio. 
Me/ Mi, a small parameter. In this linnit, B ~ a,, E ~ -(3/2)a, 
and the dispersion relation simphfies to 



1 



2A 



■-AD+C 



= 0. 



(51) 



The first factor leads to the familiar Alfven-wave dispersion 
relation: 

Lo = ±.k\\VA. (52) 

It is not hard to verify that this branch corresponds to fluctu- 
ations of (j) and A||, but not of Thus, the Alfven-wave 
dispersion relation in the fc^p, limit is unchanged from the 
MHD result. This is expected (and well known) since this wave 
involves no motions or forces parallel to the mean magnetic 
field. The wave is undamped and the plasma dispersion func- 
tion (which contains the wave-particle resonance effects) does 
not appear in this branch of the dispersion relation. To higher 
order in k±pi, however, the Alfven wave is weakly damped; 
taking the high-beta result (/3, 1), derived in Appendix D.l, 
the damping of the Alfven wave in the limit k_i_pi <C 1 is 



7 = 



(53) 



The second factor in equation (51) represents the slow-wave 
solution of the dispersion relation. This involves motions and 
forces along the magnetic field fine (perturbations of , but 
not of A II ) and, unlike in the MHD coUisional limit, is damped 
significantly (Barnes 1966; Foote and Kulsrud 1979). The 
plasma dispersion function enters through A, C, and D; to fur- 
ther simplify the expression for the slow wave, we consider the 
high- and low-beta hmits. 

In the high-beta limit, /3, 1 , the argument of the plasma dis- 
persion function for the ion terms will be small, = Tdj s/Wi ^ 
0{l/ fii) (verified by the outcome), and we can use the power- 
series expansion (Fried and Conte 1961) 



(54) 



to solve for the complex frequency analytically. The electron 
terms may be dropped because = CiiToi/Toey^'^i'ne/'niy/^ <C 
^i. Then we can approximate A ~ 1 + 7o;/7oe, C ~ i\pK^i, 
D ~ 2y'7r^, . The dispersion relation reduces to 



whose solution is ^/ = -i/ y^A', or 



UJ=-l 



\VA 



(55) 



(56) 



This frequency is purely imaginary, so the mode does not prop- 
agate and is strictly damped, in agreement with Foote and Kul- 
srud (1979). Note that ^, ^ 0(1/ confirming the a priori 
assumption used to derived this result. 

In the low beta limit, /?; ^ 1, we shall see that the phase 
velocity of the slow wave is of the order of the sound speed 
Cs = (TQe/miY/^. The electrons then move faster than the wave 
and we can drop all terms involving electron plasma dispersion 
functions because ^ Cv/v,/,^ = (m(./2m,)'/^ ^ 1. If we further 
assume that T^e ^ To,, then the ions are moving slower than the 
sound speed, so we have ~ Cs/v,;,, = {TQe/lToi)^!'^ » 1. Ex- 
panding the plasma dispersion function in this hmit gives (Fried 
and Conte 1961) 



1 

W 

Using this expansion in A ~ \+£^iZ{^i)+TQi/Toe, C ~ (,iZ(^i) and 
D ~ 2^,Z(^i), we find that the slow- wave part of equation (51) 
now reduces to A = 0, or 



(57) 



Ik 

Toe 



k\\v,h, 



+ 1\/tT' 



k\\\v,hi 



g-(w/t||V,ft()^ =0. 



(58) 



Assuming weak damping, to be checked later, we can solve for 
the real frequency and damping rate by expanding this equation 
about the real frequency. Solving for the real frequency from 
the real part of equation (58) gives 

w = ±^||c,. (59) 

This is the familiar ion acoustic wave. Solving for the damping 
gives 



_ ( i£l 1 ^-Toe/2Toi 



(60) 



This solution agrees with the standard solution for ion acous- 
tic waves (see, e.g., Krall and Trivelpiece 1986, §8.6.3) in the 
hmit k^Xj)g <C 1. Note that the a priori assumptions we made 
above are verified by this result. 
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In summary, the gyrokinetic dispersion relation in the long- 
wavelength limit, k^pi <^ 1, separates neatly into an Alfven- 
wave mode and a slow-wave mode, while the fast wave is or- 
dered out by the gyrokinetic approximation. We have seen here 
that slow waves are subject to collisionless Landau damping, 
even in the long- wavelength limit, k±pi <C 1. Therefore, if the 
scale of turbulent motions falls below the mean free path, the 
slow mode should be effectively damped out, particularly for 
high-beta plasmas. In contrast, the Alfven waves are undamped 
down to scales around the ion Larmor radius. The linear damp- 
ing of the Alfven waves at these scales is worked out in Appen- 
dices (D.l) and (D.2), where we present the high- and low-beta 
limits of the gyrokinetic dispersion relation including the ef- 
fects associated with the finite Larmor radius. The nature of the 
turbulent cascades of Alfven and slow waves at collisionless 
scales is discussed in more detail in Schekochihin et al. (2006). 

2.6.2. Short-Wavelength Limit 

At wavelengths small compared to the ion Larmor radius, 
k±Pi ^ 1, the low-frequency dynamics are those of kinetic 
Alfven waves. It is expected that, while the Alfven- wave cas- 
cade is damped around k_\_pi ~ 1, some fraction of the Alfven- 
wave energy seeps through to wavelengths smaller than the ion 
Larmor radius and is channeled into a cascade of kinetic Alfven 
waves. This cascade extends to yet smaller wavelengths until 
the electron Larmor radius is reached, k±pe ^ 1 , at which point 
the kinetic Alfven waves Landau-damp on the electrons. 

IntheUmitA:x/9, » 1, k±pe 1, wehavero(a,),ri(a,) —^ 
and Toiae) ~ Tiiae) ~ 1, whence B ~ 1, £ ~ -1. We assume a 
priori and will verify later that ~ Oik±pe) <C 1, so the elec- 
tron plasma dispersion functions may be dropped to lowest or- 
der in k±Pe. The gyrokinetic dispersion relation is then 



-A -hi 



uj' J f3i 

where A ~ 1 + Toi/Toe- The solution is 



= A, 



(61) 



(62) 



VA+2/(i + Wroi) 

This agrees with the kinetic-Alfven-wave dispersion relation 
derived in the general plasma setting (see, e.g., Kingsep et al. 
1990). Note that, for this solution, ^ 0{k±pe) as promised. 

In order to get the (small) damping decrement of these waves, 
we retain the electron plasma dispersion functions: these are 
approximated by Z(^(,) ~ ;\/7r- Then A ~ 1 + {TQj/T{)e){\ + 
iV^S.e), C ~ -/y^^e and D ~ i2(Toe/Toi)y/n£,e. Expanding the 
resulting dispersion relation around the lowest-order solution 
[equation (62)], we get 



7 = -%\\va 



k^ fl2 



1/2 



1- 



1 i+(i+ro,/ro,)A- 

2 [l + (l + 7b,/ro,-)A/2]' 



(63) 



The transition between the long-wavelength solutions of the 
previous section and the short-wavelength ones of this section 
is treated (in the anlytically tractable limits of high and low /?,) 
in Appendices (D.l) and (D.2). 



3. NUMERICAL TESTS 

' Code development continues with support from the Department of Energy's Center for Multiscale Plasma Dynamics. 



Gyrokinetic theory is a powerful tool for investigating non- 
linear, low-frequency kinetic physics. This section presents the 
results of a suite of linear tests over a wide range of the three 
relevant parameters: the ratio of the ion Larmor radius to the 
perpendicular wavelength k^p,; the ion plasma beta, or the ra- 
tio of ion thermal pressure to magnetic pressure, /3, ; and the ion 
to electron temperature ratio Toi/Toe- We compare the results of 
three numerical methods: the gyrokinetic simulation code GS2, 
the linear collisionless gyrokinetic dispersion relation, and the 
linear hot-plasma dispersion relation. 

For a wide range of parameters, we present three tests of the 
code for verification: §3.2 presents the frequency and damping 
rate of Alfven waves; §3.3 compares the ratio of ion to electron 
heating due to the linear collisionless damping of Alfven waves; 
and §3.4 examines the density fluctuations associated with the 
Alfven mode when it couples to the compressional slow wave 
around k±pi ~ 0(1). The effect of collisions on the collisionless 
damping rates is discussed in §3.5. The breakdown of gyroki- 
netic theory in the limit of weak anisotropy ~ k± and high 
frequency co ~ fi, is demonstrated and discussed in §3.6. 



3.1. Technical Details 

GS2 is a publicly available, widely used gyrokinetic simu- 
lation code, developed' to study low-frequency turbulence in 
magnetized plasmas (Kotschenreuther et al. 1995; Dorland et 
al. 2000). The basic algorithm is Eulerian; equations (25)- 
(28) are solved for the self-consistent evolution of 5D distri- 
bution functions (one for each species) and the electromag- 
netic fields on fixed spatial grids. All linear terms are treated 
implicitly, including the field equations. The nonlinear terms 
are advanced with an expUcit, second-order accurate Adams- 
Bashforth scheme. 

Since turbulent structures in gyrokinetics are highly elon- 
gated along the magnetic field, GS2 uses field-line-following 
Clebsch coordinates to resolve such structures with maximal 
efficiency, in a flux tube of modest perpendicular extent (Beer, 
Cowley, and Hammett 1995). Pseudo-spectral algorithms are 
used in the spatial directions perpendicular to the field and for 
the two velocity space coordinate grids (energy v^/2 and mag- 
netic moment v]_/2B) for high accuracy on computable 5D 
grids. The code offers wide flexibility in simulation geome- 
try, allowing for optimal representation of the complex toroidal 
geometries of interest in fusion research. For the astrophysical 
applications pursued here, we require only the simple, periodic 
slab geometry in a uniform equilibrium magnetic field with no 
mean temperature or density gradients. 

The linear calculations of collisionless wave damping and 
particle heating presented in this section employed an antenna 
driving the parallel component of the vector potential A|| (this 
drives a perpendicular perturbation of the magnetic field). The 
simulation was driven at a given frequency uja and wavenum- 
ber ka by adding an external current j^^a into the parallel Am- 
pere's law [equation (27)]. To determine the mode frequency 
uj and damping rate 7, the driving frequency was swept slowly 
(uja/iLJa <C 7) through the resonant frequency to measure the 
Lorentzian response. Fitting the curve of the Lorentzian recov- 
ers the mode frequency and damping rate. These damping rates 
were verified in decaying runs: the plasma was driven to steady 
state at the resonant frequency = cj; then the antenna was 
shut off and the decay rate of the wave energy measured. 
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Fig. 3. — The normalized real frequency ijj/k\\VA (left) and damping rate 7A|| (right) vs. fcxP; for a temperature ratio T^i/Tf^ = 100 and ion plasma beta values 
Pi = 10, 1,0.1,0.01. Plotted are numerical solutions to the gyroldnetic dispersion relation (dashed lines), numerical solutions to the hot-plasma dispersion relation 
(dotted hues), and results from the GS2 gyrokinetic code (open squares). 



The ion-to-electron heating ratio was determined by driving 
the plasma to steady state at the resonant frequency Wj, = w and 
calculating the heating of each species using diagnostics based 
on both forms of the heating equation (39). In all methods, a 
reahstic mass ratio is used assuming a hydrogenic plasma. The 
Unear GS2 runs used a single k±_ mode and 16 points in the 
parallel direction. For most runs a velocity space resolution of 
20 X 16 points was adequate; the ion and electron heating ra- 
tio runs required higher velocity space resolution to resolve the 



heating of the weakly damped species, with extreme cases re- 
quiring up to 80 X 40 points in velocity space. 

In what follows, the linear results obtained from GS2 are 
compared to two sets of analytical solutions: 

1. Given the three input parameters k_\_pi, /?,, and 7b,/roe, 
the linear, colUsionless gyroldnetic dispersion relation [equa- 
tion (41)] is solved numerically using a two-dimensional New- 
ton's method root search in the complex frequency plane, ob- 
taining the solution uj = uiGK(k±pi,pi, 7b,/7be). 
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2. The hot-plasma dispersion relation (see, e.g., Stix 1992) 
is solved numerically (also using a two-dimensional New- 
ton's method root search) for an electron and proton plasma 
characterized by an isotropic Maxwellian with no drift ve- 
locities (Quataert 1998). To obtain accurate results at high 
kA_Pi, it is necessary that the number of terms kept in the 
sums of Bessel functions appearing in the hot-plasma dis- 
persion relation is about the same as k_\_pi. The Unear hot- 
plasma dispersion relation depends on five parameters: k_\_pi, 
the ion plasma beta /?,, the ion to electron temperature ratio 
Toi/Toe, the ratio of the parallel to the perpendicular wave- 
length k\\/k±, and the ratio of the ion thermal velocity to 
the speed of light v,hjc. Hence, the solution may be ex- 
pressed as oJ =LJHp{k±pi,0i,TQi/T()e,ki\/k±,v,h,/c). The hot- 
plasma theory must reduce to gyrokinetic theory in the limit 
of < k± and Vthjc < 1, i.e., a7flp(A;_LP,-,/3i,7b,/7be,0,0) = 

iSoKikjLPi, A', Toi/Toe). 
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Fig. 5. — Ratio of the ion to electron heating Pi/Pf. Results are shown for a 
temperature ratio of Tqi/Tq^ = 100 and values of ion plasma beta ,£i, = 1.10. Plot- 
ted are derived values using the gyrokinetic dispersion relation (dashed lines) 
and the hot-plasma dispersion relation (dotted lines); results from the GS2 gy- 
rokinetic initial-value code show good agreement over nearly five orders of 
magnitude (open squares). 



3.2. Frequency and Damping Rates 

The frequency and collisionless damping rates for the three 
methods are compared for temperature ratios Toi/T^e = 100 and 
Toi/Toe = 1 and for ion plasma beta values /3, = 10, 1,0.1,0.01 
over a range of k±pi from 0.1 to 100. The temperature ratio 
Toi/Toe = 100 is motivated by accretion disk physics and the 
temperature ratio Tbi/Tbe = 1 is appropriate for studies of the 
interstellar medium and the solar wind. The real frequency 
ui and the damping rate 7 are normalized to the Alfven fre- 
quency ^||Va. The hot-plasma calculations in this section all 
have k\\/k± = 0.001 and V;a,/c = 0.001. The number of Bessel 
functions used in the sum for these results was 100, so the re- 
sults will be accurate for k±pi < 100. 

Figure 3 presents the results for the temperature ratio 



Tm/Toe = 100; Figure 4 for the temperature ratio Tqi/Tqc = 1- The 
results confirm accurate performance by GS2 over the range of 
parameters tested. 

3.3. Ion and Electron Heating 

An important goal of our nonhnear gyrokinetic simulations 
to be presented in future papers is to calculate the ratio of ion to 
electron heating in collisionless turbulence (motivated by issues 
that arise in the physics of accretion disks; see Quataert 1998; 
Quataert and Gruzinov 1999). Using the heating equation (39), 
the solutions of the hnear collisionless dispersion relation can 
be used to calculate the heat deposited into the ions and elec- 
trons for a given linear wave mode. These results for ion to 
electron power are verified against estimates of the heating from 
the hot-plasma dispersion relation and compared to numerical 
results from GS2 in Figure 5. The power deposited into each 
species Pj is calculated by GS2 using both forms of the heat- 
ing equation (39) (neglecting interspecies collisions). Here we 
have plotted the ion to electron power for a temperature ratio of 
7b(/7be = 100 and ion plasma beta values of fit = 1,10. TheGS2 
results agree well with the linear gyrokinetic and linear hot- 
plasma calculations over five orders of magnitude in the power 
ratio. 
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Fig. 6. — Electron density fluctuations for a plasma with a temperature ratio 
of Toi /roj = 1 and values of ion plasma beta /3, = 0.1, 1, and 10 using the gyroki- 
netic dispersion relation (dashed lines) and the hot-plasma dispersion relation 
(dotted lines). The fractional electron density fluctuation \Sne\/noe is normal- 
ized here by the total electron velocity fluctuations relative to the Alfven speed 
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3.4. Density Fluctuations 

Alfven waves in the MHD limit (at large scales) are incom- 
pressible, with no motion along the magnetic field and no asso- 
ciated density fluctuations. However, as Alfven waves nonlin- 
early cascade to small scales and reach k±pi ^ 1, finite-Larmor- 
radius effects give rise to non-zero parallel motions, driving 
density fluctuations. Here we compare the density fluctuations 
predicted by the linear gyrokinetic dispersion relation with that 
from hot-plasma theory. Figure 6 compares the density fluctua- 
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tions for a plasma with temperature ratio Toi/Toe = 1 and values 
of ion plasma beta /3, = . 1 , 1 , and 1 0, parameter values relevant 
to the observations of interstellar scintillation in the ISM. These 
results demonstrate that the fractional electron density fluctua- 
tions I StiQe I /noe, normalized by the electron velocity fluctuation, 
peak near the ion Larmor radius as expected. 

3.5. Collisions 

Gyrokinetic theory is valid in both the collisionless and col- 
Usional limits. To demonstrate the effect of collisionality — 
implemented in GS2 using a pitch-angle scattering operator 
on each species with a coefficient — Figure 7 presents the 
measured linear damping rate in GS2 as the collision rate 
is increased. Parameters for this demonstration are /3, = 10, 
Tm/Toe = 100, and k±pi = 1.414. The collision rates for both 
species are set to be equal, v = vu = v^e, and interspecies colli- 
sions are turned off. The figure clearly demonstrates that in the 
colUsionless Umit, v <^ijj, the damping rate due to colhsionless 
processes becomes independent of v. As the collision rate is 
increased, heating via collisionless Landau damping becomes 
less effective and the measured damping rate decreases; this is 
expected because, in the MHD limit where collisions dominate, 
Alfven waves are undamped. As discussed in §2.5, however, 
all heating is ultimately colUsional because collisions are nec- 
essary to smooth out the small-scale structure in velocity space 
produced by wave-particle interactions. A minimum collision 
rate must be specified for the determination of the heating rate 
to converge. In the low velocity-space resolution runs (10 x 8 
in velocity space) presented in Figure 7, for vl{k\\v/C) < 0.01, 
the heating rate did not converge accurately. Increasing the 
velocity-space resolution lowers the minimum threshold on the 
collision rate necessary to achieve convergence. 
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Fig. 7. — Plot of the damping rate 7 determined in linear runs of GS2 (open 
squares) vs. the coUision rate v normalized by v^. In the colhsionless limit, 
fc" ^ w, the damping rate is independent of the colhsion rate as expected. As 
the colhsion rate is increased, collisionless processes for damping the wave are 
less effective and the damping rate diminishes. 
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of applicability of the gyrokinetic solution (dashed line) are shown as the lat- 
ter deviates from the hot-plasma solution (dashed line) when Lu/Qj — > 1. Also 
plotted (solid line) is the value of kii /k± as a function of ui/Qi. 



3.6. Limits of Applicability 

The gyrokinetic theory derived here is valid as long as three 
important conditions are satisfied: (1) <C (2) a; <C O,-, and 
(3) Vih^ ^ c (the nonrelativistic assumption is not essential, but 
it is adopted in our derivation). As discussed at the end of §2.1, 
the gyrokinetic formaUsm retains the low-frequency dynamics 
of the slow and Alfven waves and collisionless dissipation via 
the Landau resonance, but orders out the higher-frequency dy- 
namics of the fast MHD wave and cyclotron resonances. 

A demonstration of the breakdown of the gyrokinetic approx- 
imation when these Umits are exceeded is provided in Figure 8 
for a plasma with 7o,/7o(. = 1, = 1, and k±pi = 0.1. Here, we 
increase the ratio of the parallel to perpendicular wavenumber 
k\\/k± from 0.001 to 100, and solve for the frequency using 
the hot-plasma dispersion relation. The frequency increases 
towards the ion cyclotron frequency as the wavenumber ratio 
approaches unity. The frequency uj/k\\VA in gyrokinetics is in- 
dependent of the value of /: | , so we compare these gyrokinetic 
values with the hot-plasma solution. Figure 8 plots the normal- 
ized real frequency uj/k^^VA and damping rate |7|/^||Va against 
the ratio of the real frequency to the ion cyclotron frequency 
Lo/fli. Also plotted is the value of k^^/k± at each ui/fli. At 
uj/r!.i ^ 0.01 and k\^/k± ^ 0.1, the damping rate deviates from 
the gyrokinetic solution as cyclotron damping becomes impor- 
tant. The real frequency departs from the gyrokinetic results at 
uj/il, ^0.1 and /k± ^ 1 . It is evident that gyrokinetic theory 
gives remarkably good results even when k\^ /k± ~ 0. 1 . 

4. CONCLUSION 

This paper is the first in a series to study the properties 
of low-frequency anisotropic turbulence in collisionless astro- 
physical plasmas using gyrokinetics. Our primary motivation 
for investigating this problem is that such turbulence appears to 
be a natural outcome of MHD turbulence as energy cascades to 



12 



HOWES, COWLEY, DORLAND, HAMMETT, QUATAERT, AND SCHEKOCHIHIN 



small scales nearly perpendicular to the direction of the local 
magnetic field (see Figure 1). Gyrokinetic turbulence may thus 
be a generic feature of turbulent astrophysical plasmas. 

Gyrokinetic s is an approximation to the full theory of col- 
lisionless and collisional plasmas. The necessary assump- 
tions are that turbulent fluctuations are anisotropic with parallel 
wavenumbers small compared to perpendicular wavenumbers, 
^11 k±, that frequencies are small compared to the ion cy- 
clotron frequency, u <C J^,-, and that fluctuations are small so 
that the typical plasma or field line displacement is of order 
Oikr^). In this limit, one can average over the Larmor motion 
of particles, simpUfying the dynamics considerably. Although 
gyrokinetics assumes uj ^ il/, it allows for k±pi ^ 1, i.e., wave- 
lengths in the direction perpendicular to the magnetic field can 
be comparable to the ion Larmor radius. On scales k±pi < 1, 
the gyrokinetic approximation orders out the fast MHD wave, 
but retains the slow wave and the Alfven wave. Gyrokinet- 
ics also orders out cyclotron-frequency effects such as the cy- 
clotron resonant heating, but retains coUisionless damping via 
the Landau resonance.^" It is worth noting that reconnection 
in the presence of a strong guide field can be described by gy- 
rokinetics, so the current sheets that develop on scales of less 
than or equal to p, in a turbulent plasma are self-consistently 
modeled in a nonlinear gyrokinetic simulation. The enormous 
value of gyrokinetics as an approximation is threefold: first, it 
considerably simplifies the linear and nonlinear equations; sec- 
ond, it removes the fast cyclotron time scales and the gyroangle 
dimension of the phase space; and third, it allows for a simple 
physical interpretation in terms of the motion of charged rings. 

In this paper, we have presented a derivation of the gyroki- 
netic equations, including, for the first time, the equations de- 
scribing particle heating and global energy conservation. The 
dispersion relation for linear coUisionless gyrokinetics is de- 
rived and its physical interpretation is discussed. At scales 
k± p, <C 1 , the familiar MHD Alfven and slow modes are recov- 
ered. As scales such that k±pi ^ 1 are approached, there is a lin- 
ear mixing of the Alfven and slow modes. This leads to effects 
such as coUisionless damping of the Alfven wave and finite- 
density fluctuations in Unear Alfven waves. We have compared 
the gyrokinetic results with those of hot-plasma kinetic theory, 
showing the robustness of the gyrokinetic approximation. We 
also used comparisons with the analytical results from both the- 
ories to verify and demonstrate the accuracy of the gyrokinetic 
simulation code GS2 in the parameter regimes of astrophysical 
interest. We note that although our tests are linear, GS2 has al- 
ready been used extensively for nonlinear turbulence problems 
in the fusion research (e.g., Dorland et al. 2000; Jenko et al. 
2000; Rogers et al. 2000; Jenko et al. 2001; Jenko and Dorland 
2001, 2002; Candy et al. 2004; Ernst et al. 2004) and a program 
of nonlinear astrophysical turbulence simulations is currently 
underway. 

In conclusion, we briefly mention some of the astrophysics 
problems that will be explored in more detail in our future work: 

1. Energy in Alfvenic turbulence is weakly damped until 
it cascades to k^p, ^ 1. Thus, gyrokinetics can be used 
to calculate the species by species heating of plasmas 
by low-frequency MHD turbulence where the dominant 
heating is due to the Landau resonance. In future work, 
we will carry out gyrokinetic turbulent heating calcula- 



tions and apply them to particle heating in solar flares, 
the solar wind, and hot radiatively inefficient accretion 
flows (see Quataert 1998; Gruzinov 1998; Quataert and 
Gruzinov 1999; Leamon et al. 1998, 1999, 2000; Cran- 
mer and van Ballegooijen 2003, 2005, for related an- 
alytical and observational results). The use of the gy- 
rokinetic formalism, which orders out high-frequency 
dynamics, such as fast MHD waves and cyclotron heat- 
ing, is justified on the assumption that the turbulent cas- 
cade remains highly anisotropic down to scales of or- 
der ion Larmor radius, so that the fluctuation frequen- 
cies are well below the cyclotron frequency even when 
k±pi ^ 1, rendering the cyclotron resonance unimpor- 
tant to the plasma heating. 

2. In situ observations of the turbulent solar wind di- 
rectly measure the power spectra of the the magnetic 
field (Goldstein et al. 1995; Leamon et al. 1998, 1999) 
and the electric field (Bale et al. 2005) down to scales 
smaller than the ion gyroradius. Thus, detailed quan- 
titative comparisons between simulated power spectra 
and the measured ones are possible. It is worth noting, 
however, that the solar wind is, in fact, sufficiently col- 
lisionless that the gyrokinetic ordering here is not en- 
tirely appropriate and the equilibrium distribution func- 
tion Fq may deviate from a Maxwellian. Significant dis- 
tortions of Fo, however, are tempered by high-frequency 
kinetic mirror and firehose instabilities that may play 
the role of colhsions in smoothing out the distribution 
function, so the solar wind plasma may not depart sig- 
nificantly from gyrokinetic behavior. It must be kept 
in mind, though, that if these instabilities were effec- 
tive in driving a cascade to higher (and thus higher 
frequency), the gyrokinetic approximation would be vi- 
olated and further analysis required to take account of 
the smaU-scale physics. 

3. In the interstellar medium of the Milky Way, the 
electron-density fluctuation power spectrum, inferred 
from observations, is consistent with Kolmogorov tur- 
bulence over 10 decades in spatial scale (see Armstrong 
et al. 1995). Observations suggest the existence of an in- 
ner scale to the density fluctuations at approximately the 
ion Larmor radius (Spangler and Gwinn 1990). These 
observations may be probing the density fluctuations as- 
sociated with Alfven waves on gyrokinetic scales (see 
§3.4) — precisely the regime that is best investigated by 
gyrokinetic simulations. 
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We are describing here the standard version of gyrokinetics. See Qin et al. (2000) for an extended theory that includes the fast wave and high frequency modes. 
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APPENDIX 

A. DERIVATION OF THE GYROKINETIC EQUATIONS 

The nonlinear gyrokinetic equation in an inhomogeneous plasma was first derived by Frieman and Chen (1982). In this Appendix, 
we derive the nonlinear gyrokinetic equation in the simplest case where the equilibrium is a homogeneous plasma in a constant 
magnetic field, i.e., VFq = and Bq = Bqx in a periodic box. We begin with the Fokker-Planck equation and Maxwell's equations. 
The gyrokinetic ordering, which makes the expansion procedure possible was explained in §2.1. The guiding center coordinates 
and Ihe key mathematical operation — ring averaging — were introduced in §2.3. In what follows, the Fokker-Planck equation is 
systematically expanded under the gyrokinetic ordering. The minus first, zeroth and first orders are solved to determine both the form 
of the equilibrium distribution function and the evolution equation for the perturbed distribution function — the gyrokinetic equation 
(the slow evolution of the equihbrium enters in the second order and is worked out in Appendix B). At each order, a ring average 
at constant guiding center position is employed to eUminate higher orders from the equation. Velocity integration of the perturbed 
distribution function yields the charge and current densities that appear in the gyrokinetic versions of Maxwell's equations. 

A. 1 . Maxwell's Equations and Potentials 
Let us start with Poisson's law. The equilibrium plasma is neutral, Qsn^s = 0. so we have 

V ■5E = A'K'Y^qMs- (Al) 

s 

The left-hand side is 0{eBov,hJ cpi) [see equation (8)], the right-hand side is OieqiUoi), so the ratio of the divergence of the electric 
field to the charge density is 0((]~^vji, /c^). Therefore, in the limit of non-relativistic ions, the perturbed charge density is zero.^^ 
This estabhshes the condition of quasineutrality: 

" : 0. (A2) 



In Faraday's law, 

cVxSE = — — , (A3) 
at 

the left-hand side is Oie^iBo), whereas the right-hand side is 0((^VLiBq). Therefore, to dominant order the electric field satisfies 
V X (5E = 0, so the largest part of the electric field is electrostatic. The inductive electric field does, however, gives an important 
contribution to the parallel electric field. The electric and magnetic fields are most conveniently written in terms of the scalar 
potential cf) and vector potential A: 

1 (9A 

5E = -V(A ^, (5B = V X A. (A4) 

c at 

We choose the Coulomb gauge V • A = 0. Thus with the gyrokinetic ordering, to O(e^), the vector potential is 

A=A||Z+A_L=A||Z + VAxz. (A5) 

Hence, the perturbed magnetic field to O(e^) is given by 

= VA|| X z- V^Az = VA|| X i + 5B\\i. (A6) 

We will use the scalars Ay and 55 y (rather than A) in subsequent development. 
Consider next the Ampere-Maxwell law: 

„ ^ 47r 1 9(5E 

Vx^ = — + - — . (A7) 

c c ot 

The left-hand side is 0{eBoQ.i/vth,), while the second term in the left-hand side (the displacement current) is Oie^BoQiVthJc^)- The 
ratio of the latter to the former is, therefore, 0(ev^^./c^), so we can drop the displacement current and use the pre-Maxwell form of 
Ampere's law: 

4-7r 

Vx 5B = -V^A = -V^A||Z-i-V(5fi|| xi=—6j. (A8) 
A.2. The Gyrokinetic Equation 

Let us start with the Fokker-Planck equation 

df. 91 q, f ISA vxB\ df, f dfA riff^.^riff^ rAQ^ 

-7- = ^+V-V/,-l-— -V0-- — -I- '^= "37 =Csrifs,fr)+Qsifs,fs), (A9) 

at ot wis V <^ c*? c J o\ \ at J ^^^^ 

where the right-hand side is the standard Fokker-Planck integro-differential collision operator (e.g., Helander & Sigmar 2002). 

Csr(fsTfr) denotes the effect of collisions of species s on (the other) species r and Css(fs,fs) denotes like-particle colUsions. To reduce 
the clutter, we suppress the species label s in this section and denote the entire collision term by C(f,f). 

Alternatively, one can say that the divergence of the electric field is small for gyrokinetic perturbations whose wavelengths are long compared to the electron Debye 
length, Xoe- 
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The distribution function is expanded in powers of e: 

f = Fo + Sf, 6f = 6fi + 6f2 + ---, (AlO) 
where 6f„ ~ 0{e"Fo). With the ordering defined by equations (8-13), the terms in the Fokker-Planck are ordered as follows 
dFo d5f , , qf ISA vx(5B vxBo\ dFo 



dt dt " m\ c dt c c J d\ 

e I e I e I 1/e 

where we label below each term its order relative to F^vthjlo- We now proceed with the formal expansion. 

A.2.1. Minus First Order, 0(l/e) 
From equation (Al 1), in velocity variables transformed from v to (v, v_\_,6), we obtain at this order 

so the equilibrium distribution function does not depend on gyrophase angle, = ^o( v, v_l , ? ) . 

A.2.2. Zeroth Order, 0(1) 

At this order, equation (All) becomes 

. vVi + f-V^+ ^) • - = aFo,Fo). (A13) 

m\ c J o\ oO 

At this stage both Fq and 5f\ are unknown. To eliminate 5f\ from this equation (and thereby isolate information about Fq), we 
multiply equation (A13) by 1 + lnFo and integrate over all space and all velocities, making use of equation (A12) and assuming that 
perturbed quantities spatially average to zero (this is exactly true in a periodic box). We find 



2v(lnFo)C(Fo,Fo) = 0. (A14) 



It is known from the proof of Boltzmann's H Theorem that this uniquely constrains Fq to be a Maxwellian: 

where the mean plasma flow is assumed to be zero. The temperature T^it) = il/2)mvjf^ associated with this Maxwelhan varies on the 
slow (heating) time scale, them ~ Oie'^lo/vth), due to conversion of the turbulent energy into heat. In Appendix B, we determine this 
heating in the second order of the gyrokinetic expansion. The density «o does not vary because the number of particles is conserved. 
In most other derivations of gyrokinetics, Fq is not determined, although it is often assumed to be a Maxwellian. 
Substituting the solution for Fq [equation (A15)] into equation (A13) and using C(Fo,Fo) = yields 

v.-VJ/,-f^^=-v.v(i^)Fo. (A16) 

This inhomogeneous equation for 6 fi supports a particular solution and a homogeneous solution. Noting the particular solution 
6fp = -{q<p/To)Fo+0{e^Fo), the first-order perturbation is written as 6f\ = -(^^/7o)Fo+/j, where the homogeneous solution h satisfies 

where we have transformed the 9 derivative at constant position r to one at constant guiding center R. Thus h is independent of the 
gyrophase angle at constant guiding center R (but not at constant position r): h = hiR,v,v±,t). Therefore, the complete solution for 
the distribution function, after taking l-q<p/To = exp{-q(j)/To) + 0(e^) and absorbing O(e^) terms into J/i, is 

q4>(r,ty 



/ = Fo(v,e^?)exp 



+ hiR,v,v±,t) + 5f2 + --- (A18) 



To 

The first term in the solution is the equilibrium distribution function corrected by the Boltzmann factor. Physically this arises from 
the rapid (compared to the evolution of 0) motion of electrons along and ions across the field lines attempting to set up a thermal 
equilibrium distribution. However the motion of the particles across the field is constrained by the gyration and the particles are not 
(entirely) free to set up thermal equiUbrium. The second term is the gyrokinetic distribution function that represents the response of 
the rings to the perturbed fields. 
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A.2.3. First Order, 0(e) 

Plugging in the form of solution given in equation (A18) and transforming into guiding-center spatial coordinates and velocity 
coordinates (v, v±,6), the Fokker-Planck equation to this order becomes 

dh dR dh q f vx(5B\ f y dh v^ dh\ fd5f2\ q / 86 v dA\ 

where 

= vuz+—[-V(p — + xz. (A20) 

dt Bo\ c dt c J 

Note that the linearized collision operator C(h,F()) + C{F(),h) involves h and Fq both of electrons and of ions. 

To eliminate 6/2 from equation (A 19), we ring average the equation over 9 at fixed guiding center R, taking advantage of the 
fact that 6/2 must be periodic in 9. The ring averaging also ehminates the third term on the left-hand side. Indeed, for an arbitrary 
function a(r), 

(-V.,. = -.((v.,<A-).V.)__.o((v.,^(-)_^.V.)_^..(,,.,.(-)j_^.-o((|)j^.O, 

(A21) 

whence (v- Vx0)r = and (v_l • (v x (5B))r = V||(v_l • (z x 5B))r = V|| (v_l • V_lA||)r = [see equation (A6)]. Thus, the ring-averaged 
equation (A19) takes the form 

''' + (^) ~-{^] =^^Fo, {All) 



dh /dR\ dh_fdh\ _ q d{x)R 



where we have defined the gyrokinetic collision operator {dh/dt)co\\ = {C(h,Fo) + C(Fo,h))R and the gyrokinetic potential x = 4'~ 
y-A/c. Keeping only first-order contrubutions in equation (A20) and substituting for SB from equation (A6), we find that the 
ring-averaged guiding center motion is given by 

where we have used the identity (v_l(5B|| )r = - (V_l(v_l • Aj,))^. Substituting equation (A23) into equation (A22), we obtain the 

gyrokinetic equation: 

9h ^ dh c r, , ,T f dh\ q9(x)R^ ..^.s 

¥ aR^^ t^^^-'^] - (a^j,,, = t^^- ^^^^^ 

where the nonhnear effects enter via the Poisson bracket, defined by 

[{x)K,h] = X • aR = -^dY-^dX- ^^^^^ 

The gyrokinetic equation (A24) describes the time evolution of h, the ring distribution function. The second term on the left-hand 
side corresponds to the ring motion along Bq, the third term to the ring motion across Bq, the fourth term to the effect of colUsions. 
The source term on the right-hand side is the ring-averaged change in the energy of the particles. A more detailed discussion of the 
physical aspects of this equation is given in §2.4. The equilibrium distribution function Fq changes only on the slow (heating) time 
scale and is thus formally fixed (stationary) with respect to the time scale of equation (A24). The evolution of Fq is calculated in 
Appendix B. 

A.3. The Gyrokinetic Form of Maxwell's Equations 

To complete the set of gyrokinetic equations, we need to determine the electromagnetic field, which is encoded by x- To determine 
the three unknown scalars (p, A|| and 5B\\ [which relates to A^ via equations (A4) and (A6)], we use the quasineutrahty condition, 
equation (A2), and Ampere's law, equation (A8), taken at (9(e) in the gyrokinetic ordering. 

A.3.1. The Quasineutrality Condition 

The charge density needed in the quasineutrality condition, equation (A2), can be determined by multiplying the distribution 
function, equation (A18), expanded to first order 0(e), by the charge q^ and integrating over velocities. Expanding the exponential in 
the Boltzmann term and dropping terms of order (9(e^) and higher gives 



E 



-^^+qj d'yhAr + - 



Tqs J V ^. 



: 0. (A26) 



Note that the velocity integral must be performed at fixed position r, because the charge must be determined at fixed position r, not 
at fixed guiding center R. Using the ring average at constant r [equation (22)], the quasineutrahty condition can be written in the 
following form 
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A. 3.2. The Parallel Ampere's Law 

The current density is calculated by multiplying the distribution function, equation (A 1 8), expanded to first order 0(e), by ^^v and 
integrating over velocities. The Boltzmann part of the current is odd with respect to V|| and vanishes upon integration. The parallel 
component of Ampere's law, equation (A8), is, therefore 

Ajr 4-7r f 

-ViA|| = —^711 = ^ — ^ J dhv^i {K),, (A28) 

where the ring average at fixed position r appears in the same fashion as in equation (A27). 

A.3 .3 . The Perpendicular Ampere 's Law 

The perpendicular component of Ampere's law, equation (A8), is derived in an analogous manner as the parallel component in 
§A.3.2: Ampere's law is crossed with z, the Boltzmann contribution vanishes upon integration over gyrophase angle 0, and a ring 
average at fixed position r is performed. The result is 

c ^ c J 

It is straightforward to show that equation (A29) is the gyrokinetic version of perpendicular pressure balance (no fast magnetosonic 
waves). Integration by parts yields: 

Vi4^=-Vi-(5Pi, (A30) 
47r 

where the perpendicular pressure tensor is 

6P± = / d\{y±y±hs}^. (A31) 

s ■' 

In order to drive steady-state (non-decaying) turbulence, we introduce additional externally driven antenna current ja to the right 
hand sides of equation (A28) and equation (A29). 

B. DERIVATION OF THE HEATING EQUATION 

While the gyrokinetic equation (A24) determines the evolution of the perturbation to the distribution function on the intermediate 
time scale, we must go to second order, Oie^uf), in the gyrokinetic ordering to obtain the slow evolution of the equilibrium distribution 
function F(). This Appendix contains two derivations of the heating equation (39): the first is more conventional, but longer (§B. 1); the 
second employs entropy conservation and is, in a sense, more intuitive and fundamental (§B.2). We also discuss energy conservation 
in driven systems and derive the power-balance equation (40) (see §B.3). 

B . 1 . Conventional Derivation of the Heating Equation 

We begin by defining the medium-time average over a period A? long compared to the fluctuation time scale but short compared 
to the heating time scale, l/a;<cA?<C 1/ e^w. 



1 



j+Aj/2 



a{t)=— / dt'ait'). (Bl) 

J I- At 12 

The equilibrium distribution function Fo is constant at times ^ Af , so Fo = Fq. 

To determine the evolution of the equilibrium density and temperature of a species s on the heating (transport) time scale, we 
consider the full (not ring averaged) Fokker-Planck equation (A9). To demonstrate that particle conservation implies that nos is a 
constant, we integrate equation (A9) over all space and velocity, divide by system volume V, and discard all terms of order O(e^) and 
higher: 

Here we have used the conservation of particles by the coUision operator, the fact that the first-order perturbations spatiaUy average to 
zero, and an integration by parts over velocity to simpUfy the result. Performing the medium-time average [equation (B 1)] eUminates 
the 6 f2s term, leaving 

-^ = 0. (B3) 
at 

Thus, for both species s, the density no.s is constant on the heating (transport) time-scale. 

The evolution of the temperature Tq^ is calculated similarly by multiplying equation (A9) by msV^/2, integrating over all space and 
velocity, and dividing by the system volume V. Using the expansion of the distribution function, integration by parts in velocity, we 
get 

3 dTn. d f d^r f m,v^ ^ f d^r /" , f d^r f , m,v^ 

T''~^^Jt 1^1 '^^^^f^^ = j ^ j d'yqsiyE)f, + J —J d\^CUfs,fr). (B4) 
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We have again assumed / d^rS f\s = (first-order perturbations spatially average to zero). The first term on the right-hand side is the 
work done on the particles by the fields, the second is the collisional energy exchange between species. Note that collisions between 
like particles do not produce a net loss of energy for a species and thus do not appear in this expression (the operator Css integrates 
to zero). At this order, the collisions between species occur only between the MaxweUian equilibria of each interacting species; the 
standard form of this collisional energy exchange (see, e.g., Helander & Sigmar 2002) is given by 

,3 /. 2 

j dh'^CMsJr) = nostynTor-Tos). (B5) 

Expressions for the interspecies collision rate t'f can be found in Helander & Sigmar (2002). The collisional energy exchange rate 
of ions on electrons, v'^, is a factor {nie/miYl^ times smaller than the ion-ion collision rate. It is, therefore, possible to sustain a 
temperature difference between ions and electrons even though the plasma is colUsional enough to make Fq MaxweUian for each 
species. 

Splitting the electric field into potentials, E = /c)dA/dt, we can manipulate the scalar potential part into a more useful 

higher order form. After some algebra and using equation (B5), we find 

where Sfs is the entire perturbed part of the distribution function. From this equation, we can check the order of the heating rate, 

l^os^ ~ Oie^ujTos). (B7) 

We see that the variation in the equilibrium quantities is, as expected, on a time scale that is slower than the time scale of the 
variations in the fluctuating quantities. Note that this ordering is consistent with all the energy in an Alfven wave cascade becoming 
heat in a single cascade time at the driving scale. 

Sphtting the perturbed distribution function Sfs into the Boltzmann and gyrokinetic parts [equation (A18)], we obtain the following 
instantaneous form of the heating equation 



3 dTi^)^ d 
2"''^^Jt 



' c/ V -^5f2s + q,(i)h ' - ' 



— / dhqs-^hs+nosi^'E'(Tor-Tos). (B8) 



V J \ 2 ' ' " J J V 2Tos 
To obtain the heating equation (39), we take the medium-time average, defined by equation (Bl), of the above equation. The average 
of the second term on left-hand side is zero; it does not contribute to the average heating. Thus, we do not require the second-order 
perturbed distribution function Sf2s to calculate the heating. The average of the first term on the right-hand side of equation (B8) is 
the desired heating term that relates the slow-time-scale evolution of the equilibrium to the solution of the the gyrokinetic equation. 
The r integral is converted into the Rj integral by noticing that / d^r J = J d^\ J d^Rs (the velocity integration on the left is at 
constant r, while on the right it is at constant Rj). Therefore, 

y IT y ^ v,.^/z.(R.) = jd^j ^q. y-x [^s-^j\ hAR.) = Jd.J -^qs^hs. (B9) 

We must now demonstrate that the heating is ultimately colUsional [the second equality in equation (39)]. To make this connection, 
we multiply the gyrokinetic equation (A24) by T^shs/F^s and integrate over space (i.e., with respect to R,) and velocity to obtain the 
following equation 

1 /,3, /,3, /f^Zo,,^ fdK\ ^ |,3, [d^,^^,^ (BIO) 



dtj J V 2Fo, W J V Fos V*/coii J J y dt 
(for reasons that will become apparent in Appendix B.2, this will be referred to as the entropy-balance equation). Combining 
equations (B8) and (BIO) using equation (B9) then gives the collisional form of the instantaneous heating equation: 

d^Y /■ 3 fnisV^ , Tqs ,j2\\ f d^Y nosql(j? 



3 dT{)s d 
xnos— 7- + -r 
2 at at 



V 2T( 



Os 



d fd\^ (h, ( ^ ) ) +no,<(%-ro.), (BID 

coll/ 



V J Fos\ \dt 



where we have used J d^Y J d^\a(Ri) = J d^Y J d^\{a)r (this manipulation is done purely for notational cleanliness: the expressions 
under the integrals are now explicitly functions of r and v, not of R,). Under medium-time averaging, the second term on the left 
hand side of equation (Bll) again vanishes and the second equaUty in equation (39) is obtained: 

h'-^-I ^ / "''K'-lf '»'^> 

The term that has averaged out does not contribute to the net (slow-time-scale) heating because it represents the sloshing of energy 
back and forth between particles and fields (on the fluctuation time scale). On the right-hand side, the colUsional term is negative 
definite for Uke-particle and pitch-angle collision operators [these are the only relevant cases: for s = i, the ion-ion colUsions dominate; 
for s = e, the dominant terms are electron-electron collisions and the pitch-angle scattering of the electrons off the ions; all other parts 
of the collision operator are subdominant by at least one factor of (me/m,)'/^]. 

The equation (B12) for the heating on the slow time scale is sign-definite, so it is, in practice, easier to average numerically 
than the instantaneous heating equation (B8): unlike in equation (B8), calculating the average heating does not require precisely 
capturing the effect of near-cancellation of the intermediate-time-scale oscillations of the instantaneous energy transfer between 
particles and waves. The net heating is always collisional, regardless of the collision rate — when collisions are small, hs develops 
smaU scales in velocity space, typically Av ~ 0{v^l^), so that the heating is independent of the colUsion rate v. As we shall see 
below, equation (B12) relates heating to the coUisional entropy production. 
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B.2. Entropy Argument to Derive the Heating Equation 

Ignoring the ion-electron collisions, whose rate is {me/rni)^/^ times smaller than that of the ion-ion collisions, Boltzmann's H 
Theorem gives the time evolution of the entropy of the ions 5, as follows 



dSi 



dt dt J V ' 



^ / dhlnfiCuifiJd- 



(B13) 



It can be easily shown that the right-hand side is non-negative (see, e.g., Lifshitz and Pitaevskii 1981) and, therefore, that entropy 
always increases. It can also be shown that the entropy increase is zero if, and only if, the distribution function is a Maxwellian. 
Expanding /; about the Maxwellian F^i, we obtain, to order O(e^), 



dt 



FoilnFoi + il+lnFoi)5f2i + 



2Foi 



fdh f , Sfu 



d5fu 
dt 



(B14) 



where we have made use of the energy conservation properties of ion-ion collisions and of the fact that Sfu spatially averages to 
zero. Evaluating the the zeroth-order (Maxwellian) part of the integral on the left-hand side and splitting dfu into the Boltzmann and 
gyrokinetic parts [equation (A 18)], we obtain the slow evolution of temperature 



3 1 dT^i d 
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(B15) 



This result is the same as equation (Bll); the heating is now manifestly expressed as the irreversible entropy production. Under 
medium-time average, the second term on the left-hand side of equation (B15) again vanishes, so the heating equation (B12) is 
recovered. 



B.3. Energy Conservation in Driven Systems 
To determine an equation for the conservation of energy in gyrokinetics, we use Poynting's theorem 



d 
dt 



d'r 



£2 c 
— + — +— <f) c/S • (E X B) = 



dh(j+ja)-E, 



(B16) 



where ja is the current in the antenna driving the system and j is the plasma current. We shall drop the surface term (the Poynting 
flux) — this is justified, e.g., in a numerical box with periodic boundary conditions. 

From §2.1, we know that \SE\^ - 0(e'^Blvf^ /c^) and \SB\^ - O(e^Bl) [equation (8)]. Thus, in the non-relativistic limit, the 
magnetic energy dominates and we may neglect the electric field energy (this is consistent with neglecting the displacement current 
in the non-relativistic ordering). We are left with 

J<5Bp 



dt 



d'r{i+ia)-E. 



Under medium-time averaging, the left-hand side vanishes and we are left with the steady-state balance: 

d\{i+ia)-E = Q. 



I- 



(B17) 



(B18) 



On the other hand, using equations (B4) and (Bll) to calculate / d^rQ • E) = J d^r'^^J d\qs(y ■ E)/„ we can convert equa- 
tion (B17) into the following instantaneous power-balance equation 
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(B19) 



The first term on the right-hand side is the nonnegative-definite collisional entropy production [see equation (B15)], the second term 
is the external energy input. The left-hand side is the time derivative of the fluctuation energy (kinetic plus magnetic): 



ThsSfi^m 



87r 



(B20) 



Note that in the large-scale limit, appropriate for MHD, we have Alfven waves, for which hg — qs (<A)r /To* — the kinetic energy 
then becomes the E x B velocity squared as it should be in MHD. 

If we medium-time average equation (B19), the left-hand side vanishes and we have the power balance between external energy 
injection and coUisional dissipation: 
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C. DERIVATION OF THE LINEAR COLLISIONLESS DISPERSION RELATION 



The dispersion relation for a linear, coUisionless gyrokinetic system is derived beginning with the hnearized coUisionless version 
of the gyrokinetic equation (25), 

dhs dhs Qs 9(v)r, , , 

^-^^iiaF = l^-^' ^^^^ 

and the field equations (A27), (A28), and (A29). First, the electromagnetic fields and the gyrokinetic distribution function are ex- 
panded in plane waves, allowing the ring averages appearing in the equations to be written as multiphcations by Bessel functions. The 
gyrokinetic equation is then solved algebraically for the distribution function. Next, this solution is substituted into equations (A27), 
(A28), and (A29), and the integration over velocity is performed using the plasma dispersion function to simplify the parallel velocity 
integrals and modified Bessel functions to express the perpendicular velocity integrals. The condition for the existence of a solution 
to the resulting set of algebraic equations is the dispersion relation. In this Appendix, the plasma species subscript s is suppressed 
when unnecessary. 

C. 1 . Solving for the Distribution Function 

First, we decompose the electromagnetic potentials into plane wave solutions of the form a(r,f) = flexp[/(k • r-wf)] (where a 
denotes <r'> or A) and the gyrokinetic distribution function into solutions of the form /!(R, v, v_l , f ) = /zexp[/(k • R-wf)]- To solve 
equation (CI) for h, we need to express (x)r in algebraic terms. Under the plane-wave decomposition, ring averages reduce to 
multiphcations by Bessel functions. For example, the terms with the scalar potential in the definition of x = </>- v • A/c, yield 



{4>{r,t))vi = cj)e 



1 

2^ 



.^_LV_L 
.-^COS^ 



(C2) 



Here we have used the definition of the zeroth-order Bessel function (Abramowitz and Stegun 1972) and the relation between the 
position and guiding center, equation (19), noting that k - [z x v/fi] = {k^_v_\_/Cl)cos9, where 9 is the angle between kj_ and the 
particle's instantaneous Larmor radius ix\ /Vt. After further algebraic manipulations of this kind, the ring averaged potential (x)r 
can be written in terms of zeroth- and first-order Bessel functions. 



(x)f 




V||A|| \ Ji{k±_v^_/Cl) mv\5B\\ 



„;(k-R-u;f) 



(C3) 



c j k±_v_\_/Q. q Bo 

where we have used the definition 6B» = i{k± x A_l) • z. Similarly, the ring average of the distribution function at constant r becomes 



{hiR,v,v±,t))r = Jo 



k±v± 



The hnearized gyrokinetic equation can now be solved for the distribution function: 

n=— -, (x)r- 

To w-A;||V|| 

Using equation (C3), this can be written in terms of the potentials as follows 
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(C4) 
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(C6) 



C.2. Performing the Integration over Velocity 



Because the solution for the distribution function equation (C6) is a product of functions of V|| and v±, the integrals over velocity 

space, / dh = J^dv\^ j[^v±dv± Jgde, in equations (A27)-(A29) can be expressed in terms of plasma dispersion functions and 
modified Bessel functions. 

Integrals over V|| , when not innmediately completed, are written in terms of the plasma dispersion function (Fried and Conte 1961) 

1 



I dx -, 



(C7) 



where ^ = |v,;, and the integral is performed over the Landau contour from -oo to +oo below the pole at x = ^ in the complex 
plane. Using this definition, we can write 
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dv 



Vth W-A:||V|| 



(C8) 



Integrations over vj_ can be written in terms of modified Bessel functions. Three such integrals arise: 
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(C9) 
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where Iq and Ii are the modified Bessel functions, a = and we have used the relation (Watson 1966), 

/•OO 1 / \ 

dxxJ„ipx)J„iqx)e-'^^ = (g) .V..^)/4«\ 



(CIO) 



C.3. Quasineutrality Condition 

Beginning with the gyrokinetic quasineutrality condition, equation (A27), we write the ring average of the distribution function 
at constant position r as a multipUcation by a Bessel function [equation (C4)] and substitute h from equation (C6) into the velocity 
integral. This gives 
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Tos 
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(Cll) 



where, in performing the velocity integrals, we have used the definitions given in the previous section. 

C.4. Parallel Ampere 's Law 
Following a similar sequence of steps for the parallel Amprere's law, equation (A28), we have 
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(C12) 



C.5. Perpendicular Ampere's Law 

It is convenient to take the divergence of the perpendicular Ampere's law, equation (A29), before processing it in the same way as 
the two other field equations. This gives 

2vi Ji {k±v±/ns) . 
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C.6. Dispersion Relation 

Before combining the three field equations derived above to produce the dispersion relation, we specify a hydrogen plasma, 
allowing us to take no; = noc and qi = -qe = e. We now divide equation (Cll) by qfnoi/Tou equation (C12) by {Airu/ck^^iqfnoi/Toi), 
and equation (C13) by {A-k / B^qinoi. Noting two manipulations, 

klk\c^ Toi _klpfk\vl 



Attlo^ qjnoi 



and 



2 Toi 
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we arrive at an algebraic Unear system of equations that can be written succinctly in matrix form as 

= 0, 



A A-B C 
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(C15) 



where lo = u>/\k\\ \va and the definitions of the coefficients A, B, C, D, E are given in equations (42-46). Setting the determinant of 
this matrix equal to zero gives the dispersion relation for Unear, coUisionless gyrokinetics [equation (41)] 



ajA 



AB + B"^ I ( ^-AD + C 



■■(AE + BCf 



(C16) 



The left-hand side of equation (CI 6) contains two factors, the first corresponding to the Alfven-wave branch and the second to the 
slow-wave branch; the right-hand side represents a finite-Larmor-radius coupUng between the Alfven and slow modes that occurs as 
k±pi approaches unity. 



ASTROPHYSICAL GYROKINETICS 



21 



D. ANALYTICAL LIMITS OF THE DISPERSION RELATION 



The linear, coUisionless dispersion relation equation (41) harbors the plasma dispersion functions Z(^s) and the integrals of the 
Bessel functions over the perpendicular velocity FoCo!.,) and Tiias). These functions can be expanded for large and small arguments, 
allowing an analytical form of the dispersion relation to be derived in these limits. The arguments in which the expansions will be 
made are 



\k\\\v,h, 
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Thus, the natural subsidiary expansion parameters in the dispersion relation are a,, /?,, nie/mi and Toi/Toe. The limit of long perpen- 
dicular wavelength, a, 1, discussed in §2.6.1, illuminates the physical meaning of each of the factors in the dispersion relation 
through a connection to the MHD Alfven and slow modes. In the short-wavelength limit, a,- ^ 1, discussed in §2.6.2, kinetic Alfven 
waves replace the MHD modes. In this Appendix, we examine the hmits of high and low /3„ while keeping a, finite. This allows us 
to connect the large- and small-wavelength asymptotics. 

D. 1 . High Beta Limit, /?, > 1 

For /?; » 1, we use the small-argument expansion of the plasma dispersion functions, Z(^s) ~ ;\/7r, because ^, = uJ/ VA ^ 1 
£,e = {me/mi)^l^{TQi/TQe)^I^S^i 1. To ensure the latter to be true, we need T^yjToc <§; {mi/me)(3i, which is not at all very restrictive. 
We can also take a,, <C 1 because me/mi<^ 1. The coefficients of the gyrokinetic dispersion relation become 
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(D3) 
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where we have dropped all terms of order 1 and higher in nie/mi. The auxiliary function Giai) introduced in the expression for D 



will be useful below. We will see that there are two interesting limits: k_LPi ~ C(/3, ^'''*), uj ~ 0(1) and ^_Lp, ~ 0{\), uj ~ 0(/3, ''^) 
(the ordering of uJ is assumed a priori and verified by the result, in the usual fashion). 

D. 1 . 1 . The limit kj_pi - 1 //3j^* 

—1 /2 

In this ordering, a, ~ ^, ~ 0(/3, ). Expanding Fo(q;,) ~ 1 -a, and Ti{ai) ~ 1 -(3/2)q;„ we find that A ~ 0(1) and B,C,D,E ~ 

—1 /2 

0(/3, ' ). The dispersion relation becomes 



-1/2. 



where B ~ a,, E ~ -(3/2)a,- and D ~ 2j(x;-\/7i7^. This is a quadratic equation for u. Its solution is 



w = -i-— 1/ — a,- ± 



16 




(D7) 



(D8) 



which agrees with the a priori ordering uj ^ 0(1). 

-1/4 



In the hmit k±pi <C /?; , we recover, as expected, the Alfven wave with weak damping [see equations (52) and (53)]: 
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For a,- > (16/9)iy7r//3;, the frequency is purely imaginary. In the intermediate asymptotic hmit /3, <C <C 1, we have 
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Fig. D9. — The real frequency (left) and damping rate (right) of the weakly damped Alfven mode derived by numerical solution of the linear, collisionless 
gyroldnetic dispersion relation (solid line) compared to the high-beta analytical limit (dashed hne). The approximate solution consists of two solutions, valid for 

both solutions are plotted to confirm that they match. 



' 0(/3( '^'*) and fcx A' ~ 0(1) [the + branch of equation (D8) and the - branch of equation (D13), respectively]. In the intermediate limit /3,. '^'^ <g ^xPi ^ 1. 



D.1.2. The limit k^pi ^ 1 

In this ordering, a, ~ 0(1), ^, ^ 0(/3,~^). Then AjBjE 0(1), C,D 0(/3,~'). The dispersion relation now is 

(Xi ( 2 
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^ D]=E\ (D12) 



Since D ~ 2ilo^-k / (iiG{ai), this is again a quadratic equation for uj. Its solution is 



A [ri(a,)- 1]2 V A- [ri(ai)- 1]2 VV A- [ri(a,)- 1]' 

In the limit k_\_pi <C 1, the two solutions are 



-1/2 

which agrees with the a priori ordering u) ~ 0(/3, ). 



(D14) 



The first solution is the damped slow mode [equation (56)] the second solution matches the weakly damped Alfven mode in the 
intermediate limit [equation (DIO)]. 

In the limit kj_pi 1, ri(a;) 0, G(ai) (ro(,/7o;)'/^(m<;/m,)^'^^, and equation (D13) reproduces the /3, 1 limit of kinetic 
Alfven waves [see equations (62) and (63)]: 

_k±pi .k]_pj jlF f Toe nie^^^^ 



in=±^-i^J-.{^—] . (D16) 
D.1.3. Summary 

Thus, at k^p; ^ Z?"'/"*, the low-frequency weakly damped Alfven waves [equation (D9)] are converted into two aperiodic modes, 
one weakly, one strongly damped [equations (DIO) and (Dll)]. At k^pi ~ 1, the weakly damped Alfven mode and the weakly 
damped slow mode [equation (D14)] are converted into two weakly damped kinetic Alfven waves [equation (D16)]. These are 
finally damped at fcxPe ~ 1- Note that the slow mode we are referring to is the weakest-damped of many modes into which the two 
MHD slow waves and the entropy mode are converted when their parallel wavelengths exceed the ion mean free path. 

The real frequency and damping rate for A = 100 and 7b,/7bg = 100 for the branch corresponding to the weakly damped Alfven 
mode are plotted in Figure D9. 



ASTROPHYSICAL GYROKINETICS 



23 




D.2. Low Beta Limit, Pi < 1 
For Pi <C 1, it turns out that A, B,C,D,E 0(1), so the gyrokinetic dispersion relation reduces to 
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We will focus on the first factor, which corresponds to Alfven modes [the long-wavelength Umit of the second factor gives the ion 
acoustic wave, see §2.6.1]. We order w~ 0(1) and consider two interesting limits: (me/OT,)(7bi/7be) <SC A ^ 1 and /3, ^nte/mi 1, 
Toi/Toe ^ 1- The solutions in these two limits are presented below. These solutions are plotted together with numerical solutions of 
the full dispersion relation in Figure DIO. 

D.2.1. The limit (OTjOTi)(7oi/7be) < /?,■ <C 1. 

In this limit, = ZZJ/v^ » 1 and = ('»e/m,)'/^(7bi77()e)'^^Ci 1 (slow ions, fast electrons). Expanding the ion and electron 
plasma dispersion functions in large and small argument, respectively, we get 
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The resulting dispersion relation is 
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where the right-hand side is small, so we can solve perturbatively for real frequency and small damping. The result is 
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where7 = Im(w)/|/:|||vA. Note that w ^ 0(1), as promised at the outset. In the Umit a,-,ae ^ 1, equation (D20) reduces to the Alfven 
wave solution, cU = ± 1 . 

D.2.2. The limit A ~ «: 1, 7b(/7be » 1- 

In this limit, ^ {mj/me)^l^ ^ 1 and, therefore, ^ {Toi/TQe)^l'^ ^ 1 (both ions and electrons are slow). Expanding all plasma 
dispersion functions in large arguments, we get 
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The resulting dispersion relation is 



aiToiae) nii 



2a7^ 



Me 



Pi-B 



cti + — Pi 



Me 



= -i(Bnf-ai)lmiA), 



where again the right-hand side is small and solving perturbatively gives 
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